REPORT  DOOlIMFNTATION  PAGE 

AD-A230  738 


Form  Approved 
OMB  No.  0704-0188 


ited  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions  searching  existing  data  sources, 
reviewing  the  collection  cf  information.  Send  comments  regarding  inis  burden  estimate  or  any  other  aspect  of 
is  burden,  to  Washington  Headquarters  Services,  Directorate  for  information  Operations  and  Repons.  1215  Jefferson 
Office  of  Management  and  Budget,  Paperwork  Reduction  Project  (070<-0i66).  Washington.  DC  20503. 


0 _ ,  wniy  (Leave  blank). 

2.  Report  Date. 

3.  Report  Type  and  Dates  Covered. 

1990 

Proceedings 

4.  Title  and  Subtitle. 

^Transient  Signal  Extraction  in  a  Multipath  Environment 


6.  Author(s). 

|r.  J.  Vaccaro,  E.  Maragakis,  and  R.  L.  Field 


5.  Funding  Numbers. 

Program  Element  No.  61153N 
Project  No.  032Q2 


Accession  No 


OHO 


DN257015 


7.  Performing  Organization  Name(s)  and  Address(es). 

Naval  Oceanographic  and  Atmospheric  Research  Laboratory 
Stennis  Space  Center,  MS  39529-5004 


8.  Performing  Organization 
Report  Number. 

PR  90:096:244 


9.  Sponsoring/Monitoring  Agency  Name(s)  and  Address(es). 

Naval  Oceanographic  and  Atmospheric  Research  Laboratory 
Stennis  Space  Center,  MS  39529-5004 


10.  Sponsoring/Monitoring  Agency 
Report  Number. 

PR  90:096:244 


11.  Supplementary  Notes. 
OCEANS 


12a.  Distribution/Availability  Statement. 

Approved  for  public  release;  distribution  is  unlimited. 


12b.  Distribution  Code. 


13.  Abstract  (Maximum  200  words). 


Ue  consider  the  problem  of  estimating  the  arrival  times  of  overlapping  ocean- acoustic  signals  from  a  received  signal 
which  consists  of  an  unknown  deterministic  signal  along  with  scaled  and  delayed  versions  due  to  multipath  propagation 
plus  additive  noise.  Our  objective  is  to  simultaneously  determine  the  transmitted  waveform  and  the  arrival  times.  The 
proposed  algorithm  obtains  approximately  maximum  likelihood  estimates  of  the  arrival  times  and  the  parameters  which 
characterize  the  unknown  signal.  Our  assumptions  are  that  the  number  of  the  different  paths  is  known  and  that  the 
signal  must  belong  to  a  parametric  class  of  signals.  We  demonstrate  the  algorithm  on  a  class  of  signals  consisting  of 
gated  sinusoids. 


$ 


DT1C 

ELECTE 
JAN  15  1991 

8 


14.  Subject  Terms. 

(U)  Transients;  (U)  Distributed  Sensors;  (U)  Coherence:  <U)  Detection; 
£U)  Classification 


15.  Number  of  Pages. 

_ 6 _ 


16.  Price  Code. 


17.  Security  Classification 
of  Report. 

Unclassified 


18.  Security  Classification 
of  This  Page. 

Unclassified 


19.  Security  Classification 
of  Abstract. 

Unclassified 


20.  Limitation  of  Abstract. 
SAR 


■.SN  7540-01-280-5500 


Standard  Form  298  (Rev  2-89) 

Presented  fry  ANSI  Std  Z39-18 
Z98-102 


Transient  Signal  Extraction  in 
a  Multipath  Environment 

It.J.  Vaccaro  and  10.  Maragakis  R.L  Field 

Department  of  Electrical  Engineering  NOAIIL,  Code  244 
University  of  Rhode  Island  Stennis  Space  Center 


Abstract 

We  consider  the  problem  of  estimating  the  arrival  times 
of  overlapping  ocean-acoustic  signals  from  a  received  sig¬ 
nal  which  consists  of  an  unknown  deterministic  signal  along 
with  scaled  and  delayed  versions  due  to  multipath  propaga¬ 
tion  plus  additive  noise.  Our  objective  is  to  simultaneously 
determine  the  transmitted  waveform  and  the  arrival  times. 
The  proposed  algorithm  obtains  approximately  maximum 
likelihood  estimates  of  the  arrival  times  and  the  parameters 
which  characterize  the  unknown  signal.  Our  assumptions 
are  that  the  number  of  the  different  paths  is  known  and 
that  the  signal  must  belong  to  a  parametric  class  of  sig¬ 
nals.  We  demonstrate  the  algorithm  on  a  class  of  signals 
consisting  of  gated  sinusoids. 

1  Introduction 

Time  delay  estimation  is  a  well  known  problem  occurring 
frequently  in  the  fields  of  sonar,  radar  and  geophysics.  In 
this  problem  the  received  waveform  consists  of  delayed  and 
sealed  replicas  of  the  transmitted  signal.  This  is  the  result 
of  different  reflections  and  attenuation  of  the  signal  in  the 
channel. 

The  received  waveform  r (f )  can  be  described  mathemat¬ 
ically  as 

M 

r(0  =  ~  f‘)  +  n(0  ,  0  <  I  <  T  (1) 

*=i 

where  s(t)  is  the  transmitted  signal,  a»  the  attenuation  fac¬ 
tor  for  path  k,  rk  the  time  delay  for  path  k,  M  the  number 
of  different  paths  and  n(t)  a  noise  component.  In  our  de¬ 
velopment  we  assume  that  the  noise  is  white  Gaussian. 

The  classical  method  for  estimating  the  times  of  arrival 
is  correlating  the  received  waveform  with  the  transmitted 
waveform.  The  peaks  in  the  correlator  output  give  the  es¬ 
timates  of  the  arrival  times.  It  can  be  shown  that  if  the 
signals  arc  separated  in  time  by  more  than  the  duration  of 
the  signal  autocorrelation  function,  the  correlator  is  equiv¬ 
alent  to  the  MLE  1 1 1.  Other  approaches  arc  given  in  |2|, 
|3|,  and  |4). _ 

®Thi»  work  w u  supported  In  part  by  Crsnl  N00014.89-K.6003  from 
the  Navsl  Ocein  »nd  Almorpherk  R»»«»rch  Laboniory  (N'OARL) 


A  completely  different  approach  has  recently  been  pro¬ 
posed  by  Kirsteins  (.*>.  Cl.  The  basic  idea  of  this  approach 
is  to  look  at  the  problem  in  the  frequency  domain.  Since 
a  delay  in  'lie  time  domain  is  equivalent  to  multiplication 
by  an  exponential  in  the  frequency  domain,  the  frequency 
domain  problem  is  one  of  fitting  weighted  complex  expo¬ 
nentials  to  the  spectrum  of  the  received  signal.  Utilizing 
an  iterative  method  of  fitting  complex  exponentials  as  in 
|7|  and  |8|,  this  approach  provides  a  way  of  estimating  the 
times  of  arrival.  In  this  algorithm  the  number  of  different 
paths  must  be  known  and  the  spectrum  of  the  source  sig¬ 
nal  must  be  nonzero.  The  requirement  that  the  number  of 
paths  must  be  known  is  not  too  restrictive  since  in  many 
eases  the  number  of  different  paths  can  be  determined  from 
the  geometry  of  the  channel. 

All  the  above  methods  require  the  source  signal  to  be 
known.  In  our  problem  we  assume  that  the  source  signal  is 
not  known.  Our  objective  is  to  simultaneously  obtain  good 
estimates  for  both  the  delays  and  the  source  signal  with  a 
minimal  amount  of  computations.  In  our  formulation  we 
assume  that  the  source  signal  belongs  to  a  parametric  class 
of  signals  which  means  that  it  can  be  completely  deter¬ 
mined  by  a  vector  of  parameters.  A  rectangular  pulse  for 
example  can  be  completely  characterized  by  its  duration, 
its  amplitude  and  its  starting  point.  By  assuming  that  the 
source  signal  belongs  to  a  certain  class  of  signals  we  have 
to  estimate  a  much  smaller  number  of  parameters  and  the 
problem  comes  much  better  defined.  In  our  development  we 
will  assume  that  the  number  of  paths  is  also  known.  Using 
those  two  assumptions  we  wilt  develop  a  method  of  obtain¬ 
ing  approximate  maximum  likelihood  estimates  of  the  time 
delays  and  the  source  signal  parameters. 

2  The  maximum  likelihood  estima¬ 
tor 

Making  the  assumption  that  the  signal  belongs  to  a  para¬ 
metric  class  of  signals,  we  can  rewrite  equation  (I)  as 

M 

r(t)  =  X>*a(‘  -  r»;*J  4  n(r)  .  0<f<r  (2) 

k=  I 

where  0  is  the  vector  of  parameters  which  characterize  the 
source  signal.  In  the  case  that  n(f)  is  white  Gaussian  noise 
the  least  squares  estimator  is  also  the  maximum  likelihood 
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estimator  (llelstrom,  pp  199  |9|),  Therefore  the  MLE  is 
given  liy 

rT 

min  /  r[t)  Y  „,*(/  rk,0)'7<fl  (3) 

"  ' 4  .  •* »  Jit  *T'| 

We  remark  that  one  advantage  of  the  maximum  like¬ 
lihood  formulation  of  the  problem  is  that  the  ease  when 
the  signal  is  unknown  is  conceptually  the  same  as  when 
the  signal  is  known.  If  the  signal  is  known,  then  the  pa¬ 
rameter  vector  0  in  (3)  is  fixed  and  is  not  included  in  the 
minimization.  If  the  signal  is  not  known,  then  the  (3)  must 
he  minimized  with  respect  to  all  the  parameters. 

Assuming  the  integrand  in  (3)  is  zero  outside  |0,  T|,  we 
can  extend  the  limits  to  infinity.  Using  Parseval’s  theorem 
we  get  the  equivalent  expression 

min  [  -  S(w; 0)  Y"  ake~’’tu\'1du  (4) 

Tk  .  ak  J -oo  t  . 

where  /f(w)  and  S(ui;0)  are  the  Fourier  transforms  of  r(t) 
and  s(t\0),  respectively.  It  should  be  noted  that  this  ex¬ 
pression  for  the  MLE  is  only  valid  for  white  Gaussian  noise. 
However  even  if  the  noise  is  not  Gaussian  the  estimator  may 
still  be  useful  as  the  least  squares  estimator. 

The  problem  now  is  to  approximate  R(w)  by  a  weighted 
sum  of  complex  exponentials.  If  we  sample  the  frequency 
functions  with  spacing  Aw,  we  have 


plitudes  is  a  well-known  problem,  and  accounting  for 
the  weights  is  simple.  However,  constraining  the  am¬ 
plitudes  to  be  real  when  the  data  is  complex  has  ap¬ 
parently  not  been  considered  before.  We  show  in  the 
next  section  how  to  include  this  constraint. 

3  The  known  signal  algorithm 

Using  the  notation  developed  in  the  previous  section,  we 
now  consider  the  case  when  the  source  signal  is  known  and 
is  narrowband.  We  show  in  the  next  section  how  the  known- 
signal  algorithm  can  be  used  iteratively  in  the  case  when 
the  signal  is  not  known. 

When  the  signal  is  narrowband,  most  of  the  energy  of 
the  signal  is  concentrated  in  the  passband.  For  example  in 
a  gated  sinusoid  most  of  the  energy  of  the  signal  is  concen¬ 
trated  in  the  main  lobe  around  the  center  frequency.  In  this 
case,  we  do  not  have  to  include  all  frequency  points  in  the 
minimization;  we  only  have  to  include  those  which  contain 
some  signal  energy.  If  we  take  N  points  starting  at  q  corre¬ 
sponding  to  positive  frequencies  where  the  spectrum  of  the 
transmitted  signal  is  nonzero,  we  can  define  the  following 
error  function 

£|(A*.  «*)  =  "£  l*("  +  9)  ~  w(n  +  q)  £  (6) 

o=0  4=1 


where 


Ak  =  -rkAw. 


After  sampling  the  integrals  arc  approximated  by  sums  and 
the  MLE  is  given  approximately  by 


,  ™in0(  \R(n )  -  W(M)  £  ak«;Mn> 


where  ut(m;0)  =  S(mAt»/;0),  and  L  is  the  total  number 
of  points  in  the  discrete  Fourier  transform  applied  to  the 
sampled  data. 


Remarks 

Before  giving  an  algorithm  for  estimating  parameters  from 
frequency  domain  data,  we  first  give  some  features  of  the 
frequency  domain  formulation  of  the  MLE. 


Note  that  the  above  equation  is  not  equivalent  to  the  MLE 
expression  shown  in  (5)  because  it  only  corresponds  to  the 
positive  frequency  portion  of  the  spectrum.  The  conjugate 
symmetric  portion  of  the  spectrum  corresponding  to  nega¬ 
tive  frequencies  must  also  be  included  to  obtain  the  MLE 
error  expression.  We  first  define  some  notation. 

Let 

r  =  !*(*)  *(9  +  1)  *(<?  +  *- l))r 

o  =  (a,  a,  •••  aM |r 

W  =  diag{w(g)  w(q  +  l)  •  •  •  w(q  +  N  -  1)} 
e/*n  ... 

e/M*+i)  M  ...  ,./>«(«+>) 

d(A)  .  ;  :  . 

e;A,U  +  /V- 1)  e;i,(«  +  /V- 1)  ...  e/M«  +  W-l) 

P(  A)  =  WA{  A) 

Then 

£,(a,A)  =  ||r  -  P(A)a||’.  (7) 


1.  Note  that  the  delay  parameters  are  estimated  as  real 
numbers  even  when  using  sampled  data.  There  is  no 
need  for  interpolation  as  there  would  be  in  a  time- 
domain  formulation. 

2.  The  frequency-domain  formulation  is  equivalent  to 
modeling  the  spectrum  of  the  received  signal  as  a 
weighted  sum  of  complex  exponentials  with  real-valued 
coefficients.  The  complex  exponentials  do  not  ojeur 
with  conjugate  symmetry,  in  general,  The  fitting  of 
unweighted  complex  exponentials  with  complex  am- 
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This  is  the  error  expression  considered  in  |6|.  However, 
the  MLE  error  expression  must  also  include  the  conjugate 
symmetric  portion  of  the  spectrum  as  shown  below 

£(a,A)  =  ||r  -  P(A)a||*  +  ||r*  -  r(AH|*,  (8) 

where  the  superscript  *  refers  to  complex  conjugation.  Note 
that  the  vector  a  is  not  conjugated  in  the  second  term  of 
the  above  equation  since  it  is  assumed  to  be  a  real  number. 
Thus  the  formulation  In  (8)  is  equivalent  to  the  constraint 
that  the  amp'itudcs  be  real  valued. 

. awaatty  COdM 

A *ail  and/or* 
list  Sptolal 

IV'I 


l/v 


□□ 


For  any  fixed  A,  the  coefficients  a*  which  minimize  K, 
e  given  by 

<i,  l>'(\)r  -.(P"l>)  'l>"r.  (<)) 


>te  that  the  resulting  vector  a,  will  be  complex  in  general 
ibstituting  (9)  into  (7)  yields 


_ 

E,(A)  -  ||  (/  -  P(P"P)-'P")r |!5  =  ||/u(A)r||2.  (10) 


Similar  expressions  can  be  written  for  the  true  MLE 
'ror  expression  as  follows 


II 

r  j 

P(A) 

11 

r*  J 

r-(\) 

or  any  fixed  A,  the  coefficients  a  which  minimize  E  are 
iven  by 


a 


(12) 


ami  !<(j )  W  '(j)R(]). 

Note  that  5(z)  =  60  f  b,z  -(•••(  bpzp  is  a  polynomial 
whose  roots  arc  e'\  {  -  1 ,  ,  |„  order  for  the  roots  to 

be  on  the  unit  circle  we  must  impose  the  conjugate  symme¬ 
try  constraint  l>\  -  b'M  t  as  in  |f>]  and  jH|.  We  can  enforce 
the  conjugate  symmetry  condition  on  b  by  writing  b  -  Cx 
where  r  is  a  real-valued  vector  whose  elements  consist  of 
the  real  and  imaginary  components  of  b0  •  •  •  bM/7,  and  C  has 
the  form 


i 

J 

0 

0 

0 

0 

0 

0 

1 

) 

0 

0 

1 

J 

i 

- ; 

0 

0 

.  .  0 

0 

0 

0 

1 

J 

0 

0 

0 

0 

0 

0 

■  1 

"7  . 

Using  b  -  Cz,  we  can  write  the  error  as  a  function  of 
the  vector  z  as 


E,(x)  =  ||  B(BHB)  1 Y  Cz||2.  (18) 


a  =  |Rcal(P/'P)|-|Rcal(P"r)  (13) 

vherc  the  superscript  H  stands  for  complex  conjugate  trans- 
losc,  and  Rcal(  )  stands  for  the  real  part  of  a  complex  num- 
>cr.  Note  that  the  above  expression  always  results  in  real 
/alucs  for  the  amplitudes. 

If  we  define 

P(A) 

/-(A) 

ind  substitute  (13)  and  (M)  into  (11)  we  get 


G  = 


(14) 


E( A)  =  !!(/  -  G(GhG) 


-M 


=  IIGMA) 


(15) 

The  error  is  now  only  a  function  of  the  unknown  time 
delays.  Note  that  GJ-(A)  is  the  matrix  which  projects  onto 
the  orthogonal  complement  of  the  column-space  of  G. 

The  expressions  for  E\  (which  uses  half  the  spectrum) 
and  E  (which  uses  the  complete  spectrum)  are  of  the  same 
form  as  shown  in  (10)  and  (15).  It  is  shown  in  [6]  that  E\ 
can  be  written  as  a  function  of  a  new  parameter  vector  b 
as  follows 

£,(&)  =  ||£(B"JS)-,y6!i2 

where 


B" 


bo 


bo 

0 


bM 


W 


-1 


b  = 


(10) 

feo 

bM 


bo,---, bp  arc  chosen  to  satisfy 

bo  +  &i («'*')  +  •••  +bN(eiy‘M)  =0  ,  t 


!,•••.  M, 


k(q  +  M  +  1) 
Tl(q  +  M  +  2) 


fl(q+  1) 


h(q±N  -  1)  •••  K{q  +  N  -  M) 


Before  developing  an  algorithm  to  minimize  this  error  ex¬ 
pression,  we  first  look  at  E i  and  E  for  a  one-path  example 
to  gain  some  insight  into  the  relationship  between  these  two 
error  expressions. 

3.1  A  One-Path  Example  to  Compare  E 
and  E, 

The  errors  Ex  and  E  can  be  written  as  functions  of  the 
unknown  time  delays  only  as  shown  in  (10)  and  (15).  In 
thi.  section,  we  give  an  example  of  a  signal  containing  a 
single  time  delay  to  compare  E ,  and  E. 

The  transmitted  signal  in  this  example  consists  of  a  244 
Hz  gated  sinusoid  whose  duration  is  40  ms.  The  received 
signal  was  obtained  by  delaying  the  transmitted  signal  by 
50  ms.  Thus  the  error  surfaces  should  have  minima  at 
t=0.05  secs. 

In  Fig.  1,  we  show  the  received  signal  with  a  moderate 
amount  of  noise  added.  In  Fig.  2,  we  show  the  error  sur¬ 
faces  for  E  and  E \  corresponding  to  the  received  data  in 
Fig.  3.  Note  that  both  E  and  E\  have  minima  at  t  =  0.05, 
but  E j  is  a  smooth  unimodal  function,  while  E  is  a  modu¬ 
lated  sinusoid.  Our  attempts  to  minimize  E  converged  to 
some  local  minimum  unless  the  algorithm  was  initialized 
very  close  to  the  global  minimum.  On  the  other  hand,  our 
algorithm  converged  to  the  global  minimum  of  E i  for  a  wide 
range  of  initializations. 

By  looking  only  at  Fig.  2,  we  might  conclude  that  min¬ 
imizing  E\  always  gives  the  same  results  as  minimizing  the 
true  MLE  expression  E.  However,  this  is  not  the  case.  Fig. 
3  shows  the  same  received  signal  with  a  large  amount  of  ad¬ 
ditive  noise,  and  Fig.  4  shows  the  corresponding  error  sur¬ 
faces  E i  and  E.  The  error  surface  E  still  has  a  minimum  at 
t  =  0.05,  but  the  minimum  of  E\  occurs  at  l  =  0.048.  Thus 
it  seems  that  minimizing  E\  will  result  in  biased  estimates 
of  the  time  delays  at  low  signal-to-nol”*  ratio.  Neverthe¬ 
less,  we  choose  to  work  with  E\  instead  ri  E  because  it  is 
easier  to  find  the  global  minimum  of  E\  than  it  is  for  E. 
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As  shown  in  Ihc  above  example,  Ei  will  work  well  provided 
the  signal-to-noisc  ratio  is  large  enough.  The  development 
of  an  efficient,  algorithm  to  find  the  global  minimum  of  /•,’ 
for  a  wide  range  of  initial  estimates  is  an  important  open 
problem 

3.2  A  Perturbation  Expansion  Approach 
to  Minimizing  E\ 

An  iterative  algorithm  for  minimizing  E\  has  been  derived 
in  lr>!-|8j,! lOj.  However,  the  algorithm  described  in  these 
references  failed  to  converge  for  the  example  given  in  t he 
next  s<-<  t ion  of  this  paper.  Here  we  briefly  describe  a  new 
approach  to  minimizing  E\. 

We  begin  by  considering  the  value  of  the  error  function 
at  an  increment  Ax  away  from  a  nominal  value  x  of  the 
coefficients 

E[x  «  At.)  j|/f(z  )■  Ar)(Zt"(z*  Az)Z?(z  -t  Ax)| YCx  |J. 

(10) 

After  some  calculation,  a  first-order  perturbation  expansion 
for  E(x  i  Ax)  can  be  obtained.  The  result  is 

E(x  t  Ax)  i  \\PDy  +  QD&I)"P£y  +  P^BQMByf  (20) 

where  A  means  “equal  to  first-order  in  Ax,"  Qg  = 

~  fl)~'  I)" ,  P£  =  I  ~  PB,  and  D  is  evaluated  at 

the  nominal  parameter  vector  x.  Finally,  for  any  vectors 
U|  and  V],  we  derive  first  order  expressions  relating  A/f 
and  Ax  as  follows  (that  is,  we  derive  expressions  for  the 
matrices  and  M7) 

A/?u,=A/|Az  .  , 

A5"v,=AY,Ax. 

Substituting  the  above  expressions  into  (20)  with  u,  - 
QBV,  vi  -  PfiV-  ar,d  b=Cx  yields 

E(  x  +  Ax)  =  ))u  +  MAxj|l  (22) 

where  u  -  PBy  and  M  =  M i  +  Mj.  We  then  solve  for 
Ax  which  minimizes  the  above  expression  subject  to  the 
constraint  that  Ax  is  real  valued.  The  result  is 

Ax  =  (Rcal(M"A/)|-'Real(A/"u).  (23) 

If  x  corresponds  to  a  given  value  of  the  parameter  vec¬ 
tor  6,  then  we  replace  x  by  x  +  Ax  and  repeat  the  above 
calculations  (solve  for  a  new  Ax).  The  process  continues 
until  convergence. 

4  Simultaneous  signal  extraction 

Wc  now  turn  to  the  question  when  the  source  signal  is  also 
unknown.  In  this  case  we  also  want  to  estimate  the  vector  0. 
In  our  analysis  wc  will  assume  that  the  signal  belongs  to  the 
parametric  class  of  signals  of  gated  sinusoids  of  unknown 
frequency  and  duration.  This  class  of  signals  was  chosen 
because  of  availability  of  experimental  data  of  this  class, 
First  wc  observe  that  those  can  be  completely  described 
by  specifying  the  frequency,  the  duration  and  the  starting 


point,  The  phase  of  the  sinusoid  at  the  starting  point  may 
also  be  included  with  the  above  parameters.  However  if 
the  frequency  and  the  duration  are  such  so  that  there  are 
several  cycles  of  the  sinusoid  ill  the  source  signal,  the  phase 
will  not  be  an  important  parameter.  In  order  to  reduce 
the  calculations,  this  assumption  was  made  in  our  model 
and  the  phase  was  taken  to  be  zero.  Also  note  that  the 
time  delays  are  calculated  relative  to  the  starting  point. 
Changing  the  starting  point  of  the  source  signal  will  change 
the  time  delays  by  the  same  amount.  In  our  algorithm  we 
assume  that  the  source  signal  starts  at  zero  and  wc  calculate 
the  time  delays  relative  to  that  point. 

By  making  the  above  assumptions,  the  source  signal 
extraction  has  become  the  estimation  of  the  frequency  and 
the  duration.  The  maximum  likelihood  estimator  for  the 
frequency,  duration  and  the  delays  is  equivalent  to 

rnin  |R(n)  -  w(n,f,d)  £  (24) 

n=0  t=i 

where  /  and  d  arc  the  frequency  and  the  duration  respec¬ 
tively.  This  expression  is  not  easy  to  minimize.  It  is  much 
more  complicated  than  a  fitting  of  exponentials  since  te  is 
also  nonlinear  in  f  and  d.  A  multidimensional  search  is 

practically  impossible  because  of  its  computational  inten¬ 
sity. 

Our  algorithm  reduces  the  computations  by  breaking 
the  problem  down  and  solving  for  the  parameters  and  the 
delays  as  follows: 

1.  Obtain  initial  estimates  f°,  d°  of  the  unknown  fre¬ 
quency  and  duration. 

2.  Use  /'  and  d'  in  the  known-signal  algorithm  to  esti¬ 
mate  the  time  delays  rj, 

3.  Using  the  estimated  time  delays  r^,  calculate  new  val- 
ues  /'+,i  for  the  frequency  and  duration. 

4.  Check  for  convergence;  return  to  2. 

The  questions  that  we  now  have  to  address  are  how  to 
estimate  /'  and  and  how  to  obtain  the  first  estimates  of 

f°  and  <P. 

Consider  the  error  function  for  a  given  set  of  time  delays: 

~  52  \P(n+<i)-w(n+q-,/,d)52oke*Xlin*')\i  (25) 

k»l 

E(f is  a  nonlinear  function  of  /  and  d  making  the 
minimization  very  hard  analytically.  Notice  that  the  du¬ 
ration  d  is  a  discrete  variable  since  the  signal  is  sampled. 
Because  of  this,  a  gradient  based  algorithm  of  minimizing 
E(f,d)  would  not  work  for  d.  The  advantage  of  this  is  that 
only  a  finite  number  of  values  for  d  is  possible.  E(/,d)  is 
minimized  as  follows: 

1 .  Do  a  search  over  values  of  d  near  the  previous  estimate 

of  d. 
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2.  For  each  value  of  d  find  the  best  value  of  /  using  a 
gradient  based  technique. 

3.  Repeat  2  until  a  minimum  is  found. 

•  ote  that  during  tin?  minimization  E(f,d)  will  have  to  he 
valualed  several  times.  Each  evaluation  is  made  faster  by 
he  fact  that  u ;(•;/, d)  can  be  computed  by  an  FFT. 

The  last  question  is  how  to  obtain  the  initial  estimates 
0  and  d°.  Those  estimates  do  not  have  to  be  very  accurate, 
n  fact  they  can  be  quite  crude  since  they  are  only  used  to 
nitializc  the  algorithm. 

For  our  experimental  arrangement  the  geometry  of  the 
hannel  gives  us  that  the  first  two  paths  will  be  the  least 
ittenuatcd  and  will  probably  be  overlapping.  After  filtering 
he  reflections  from  the  first  two  paths  will  lie  over  the  noise. 
Taking  advantage  of  this  we  can  gel  estimates  of  /  and  d 
us  follows. 

•  Obtain  the  envelope  of  the  signal  as  tollows:  t or  earn 
point  take  the  maximum  amplitude  over  the  next  20 
points  (generally  choose  a  number  of  points  sufficient 
to  cover  a  period  of  the  sinusoid). 

•  Take  the  initial  estimate  of  the  duration  to  be  (Q-L- 

20)/2. 

•  The  initial  estimate  of  the  frequency  is  found  using  a 
standard  frequency  estimation  algorithm  on  the  data 
points  between  L  20  and  Q. 

This  completes  the  description  of  our  algorithm.  In  the 
next  section  wc  will  demonstrate  the  algorithm  on  both  real 
and  simulated  data. 

5  Example  with  Experimental  Data 

In  our  experiments  the  geometry  of  the  channel  gives  us 
four  different  paths  with  the  first  two  less  attenuated  than 
the  other  two.  This  situation  arose  in  ocean  acoustic  signals 
reflected  on  the  surface  and  the  bottom  of  the  ocean.  When 
the  data  was  collected,  a  hydrophone  near  the  transmitter 
recorded  the  actual  source  signal.  The  transmitted  signal 
was  a  gated  sinusoid  of  frequency  244  Hz  and  duration  40 
ms. 

A  record  of  received  data  is  shown  if  Fig.  5.  Using 
the  known-signal  algorithm  presented  in  section  3  of  this 
paper,  wc  estimated  the  time  delays  and  amplitudes  of  the 
four  paths.  Wc  then  constructed  an  estimate  of  the  received 
signal  using  our  estimates  of  the  delays  and  amplitudes  as 
well  as  the  known  source  signal.  The  actual  received  signal 
and  our  reconstructed  signal  our  shown  together  in  Fig.  6, 
which  shows  that  the  estimates  provide  a  good  fit  to  the 
data. 

Next,  wc  applied  the  unknown  signal  algorithm  pre¬ 
sented  in  the  previous  section  to  try  to  estimate  the  pa¬ 
rameters  of  the  transmitted  signal  as  well  as  the  delays  and 
amplitudes  in  the  received  signal.  The  estimates  of  both  the 
known  signal  and  unknown  signal  algorithms  are  shown  in 
the  Table  1 . 


Wc  believe  that  the  estimates  obtained  by  both  algo¬ 
rithms  would  be  improved  by  minimizing  error  surface"/-,' 
instead  of  /-,',  We  also  remark  that  the  estimates  of  the 
amplitudes  are  extremely  sensitive  to  the  values  of  the  cs(i 
mated  delays,  however  the  amplitudes  arc  always  rah  ulaied 
to  give  a  least-squares  fit  to  the  data  for  a  given  set  of  time 
delays. 


True  parameter 

K.S.  estimate 

U.S.  estimate 

Frequency 

244.0 

- 

239.3 

Duration 

0.0400 

- 

0.0365 

ri 

- 

6.1550 

0.1550 

U 

- 

0.1640 

0.1657 

- 

0.1920 

0.1902 

U 

- 

0.1995 

0.2004 

U| 

- 

0.6448 

-0.2657 

- 

-0.6463 

-0.7162 

“3 

- 

-0.4237 

0.0304 

- 

-0.0207 

0.4871 

Table  1:  Estimates  of  delays  and  amplitudes  using  the 
known  signal  (K.S.)  and  unknown  signal  (U  S.)  algorithms. 
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Figure  1  A  received  signal  consisting  of  a  single  path  with 
a  moderate  amount  of  additive  noise 
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Figure  2  The  error  surfaces  E  and  E|  corresponding  to  the 
received  signal  in  Fig-  1-  The  smooth  curve  is  E\  and  the 
oscillating  curve  is  E. 
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Figure  .1  A  received  signal  consisting  of  a  single  path  with 
a  large  amount  of  additive  noise 
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Figure  4:  The  error  surfaces  E  and  E,  corresponding  to  the 
received  signal  in  Fig.  2  The  smooth  curve  is  Ej  and  the 
oscillating  curve  is  E. 


Figure  o.  A  record  of  experimental  data  with  four  overlap¬ 
ping  paths  moderate  amount  of  additive  noise. 


Figure  6:  The  received  signal  of  Fig  5  together  with  the  re¬ 
constructed  signal  obtained  by  the  known  signal  algorithm. 
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